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ABSTRACT 

Astrophysical plasmas are subject to a tight connection between magnetic fields and the diffusion of particles, which 
leads to an anisotropic transport of energy. Under the fluid assumption, this effect can be reduced to an advection- 
diffusion equation, thereby augmenting the equations of magnetohydrodynamics. We introduce a new method for solving 
the anisotropic diffusion equation using an implicit finite-volume method with adaptive mesh refinement and adaptive 
time-stepping in the ramses code. We apply this numerical solver to the diffusion of cosmic ray energy and diffusion 
of heat carried by electrons, which couple to the ion temperature. We test this new implementation against several 
numerical experiments and apply it to a simple supernova explosion with a uniform magnetic field. 
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1. Introduction 

Diffusive processes are relevant at various scales in astro- 
physical systems from the diffu sion of cosmic rays (CRs) 
in the interstellar medium (e.g. IStroner et al.ll20d^ to the 
the diffusion of heat and CRs in c l uster of galaxies (e.g. 
Rosn^er fc TuckeilFl989l : lFabianl[l994lNaravan fc MedvedevI 

200lh . As the diffusion of these charged particles oper¬ 
ate in a magnetised plasma, the diffusion process becomes 
anisotropic since it is conducted along magnetic field lines. 
Furthermore, for a fully ionised plasma, like the lighter pop¬ 
ulation of particles, electrons conduct heat much faster than 
protons, which are virtually in the limit of non-diffusion. In 
some situations, it leads to a non-thermal equilibrium be¬ 
tween ions and electrons if the collisional coupling time of 
the two temperatures is longer than the diffusion timescale 
of electrons. Thus, there could be fluctuations in electron 
temperature without the corresponding fluctuation in the 
ion temperature. 

Those processes have mostly been ignored from 
the numerical astrophysics perspective with a few ex¬ 
ceptions. In the intra-cluster medium, it has been 
shown that ions an d electrons a r e not p erfectly at ther- 
mal 


_ , uilibrium ... 

2004 iRudd fc Naeail 


(IChieze et al.l 

1 20091: 


l2004f) . Thermal and magnetohydrodynamical instabilities 
form because of the anisotropic nature of the conduc- 
tion, the so-called magne to-thermal instability (jBalbusI 


I2OOOI : iParrish et al.l 2008f) an d the heating-buoyancy 
insta- b ilitv llParr ish fc Qi iataertl 120081 : iBogdanovic et al.l 
I2OO9I : iParrish et al.l 2009(1 . and their variations with the 


adjun ction of the diffusion of CRs (|Rasera fc ChandranI 
I2OO8II . Conduction may also play an important role 
in redistributing the mechanic al energy deposited by 
activ e galactic nuclei je t s (iRuszko wski fc Begelm^ 


activ e galactic nucfei le t s IIKuszkowski &: fje gefmanl 
200^_l&ighenti^^.Mathe^ 1200,4 iBriiggen et al.l 120051 : 


McNamara fc Niilseiil l2007ll . helping to solve the cooling 


catastrophe in galaxy clusters. A similar heating source in 
the intra-cluster medium is that of comic rays deposited 
in active galactic nuclei or a t shocks and redist ributed by 
diffusion (iMiniati et al.l[200ll : lEnfilin et al.ll2011ll . 


since CR energy in the inter- 
is at equi partition with kinetic , 


1994 ICourtv Y Alimil 

_ _ - - - IWong fc SarazinI 120091 : 

Gaspari fc Chiirazovl 201^ and that the transport of 
heat produces an electron pre cursor at the virial radius 
of massive dark matter halos (|Tevssier et al.l 1199811 . The 
transport of heat eventually reduces the cooling flow 
at the core of galaxy clusters, transporting the cosmic 
shock-drive n internal energy of the ou t skirts towards 
the centre (jRuszkowski fc BegelmanI [20021 : l.liibelgas et al.l 


On galactic scales, 
stellar medium (ISM) 
thermal, and magnetic energies ( Beck fc Kraiisd 1200, 5ll . 
they may play a crucial role in the self-regulation of 
the ISM dynamics and, thus, of star formation. For in¬ 
stance, the injection of CRs within re mnants of supernovae 


lead to stronger g a .lactic-scale w i nds (l.lubelgas et ^ 
Uhhg_et_^ 2014 Hanasz etld] 12014 iBooth et alT 


2008 


2013 


Salem fc BrvanI 1201411 an d to amplihed gala ctic dynamos 


with anisotropic diffusion (iHanasz et al.ll2bb4l . Anisotropic 
conduction in the ISM affects the shape and size of 


cold clouds (iKov ama fc Inutsukall2004 iPiontek fc Ostrikeil 

|2004 IChoi &: Stonell2 012ll. and~tlie expansion of sup ernova 


I^UU4t lUhoi &: btonell/i Ul^). and the expansion 01 sup ernova 
remnants ( Tillev et al.l 1200^ iBalsara et al.l l2008af) . Also, 
CRs can penetrate deep inside cold dense cores and, as 
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a conse quence, can regulate the ionisation rates of the 
gas fe.g. ISoitzer fc TomaskolflO^ iPadovani et aLll2009li . 

One problem with the implementation of diffusion in nu¬ 
merical simulation is that the stability criterion is Atdiff = 
Aa;^/(2Ddiff), where Aa; is the cell size and Ddiff the 
diffusion coefficient. For comparison, the hydrodynami- 
cal Courant Friedrich Levy (CFL) condition is Ath = 
CAi/ (m-|-Cs), where u is the gas velocity, Cs the ffuid sound 
speed, and C < 0.8 is the Courant factor. Since the diffusion 
stability time step does not scale linearly with resolution, 
unlike the hydrodynamical time step, it could be a bottle¬ 
neck to employ an explicit diffusion scheme that must verify 
this time step condition at any time, especially for multi¬ 
scale problems where gravity shapes strong contrasts in gas 
densities and triggers refinements in resolution. 


Implicit numerical solvers in these situations could 
be favoured. They do not need to fulfil the time step 
constraint on diffusion, and any time step can be cho¬ 
sen (such as the CFL time step), at the expense of 
numerical complexity and intensity compared to an ex¬ 
plicit solver. There is a variety of numerical implemen¬ 
tations of anisotropic diffusion f rom centred symm etric 
to centred asymmetric schemes (jCfinter et al.l l2005^ us¬ 
ing various slo pe limiters to preserve th e monotonicity 
of the system (|Sharma fc Hammett! l2007f) . Amongst all 
these implementations, implicit and expl icit solvers have 
been d eveloped for astr o physi cal co des: iHanasz fc Lesch 
(2003), iParrish fc Stond l|2005^ , and iRasera fc Chandran 
( 2003 ) implemented an explicit meth o d for anisotropic 
conduction and diffusion. Balsara et al. (l2008bl) deve loped 
a semi-impl i cit so lver. lYokovama fc Shibatal (l200lll and 
iMever et al.l l|2012l) a fully implicit solver. 


_ ^ this paper, which uses the RAMSES code (iTevssietl 

I2002D . we present the first implementation of an implicit 
solver for anisotropic diffusion or conduction on adaptive 
mesh refinement (AMR) with adaptive time-stepping. This 
new numerical implementation is augmented by modelling 
a multi-temperature and energy component with tempera¬ 
ture coupling of ions and electrons. In section[51 we present 
the new numerical solver for diffusion or conduction and 
temperature coupling. In section[31 we test our method with 
several numerical experiments. We discuss future develop¬ 
ments and applications of this method in section 0] 


2. Numerical set-up 

2.1. Magnetohydrodynamics with heat and cosmic ray 
diffusion 


The set of differential equations to be solved for magneto¬ 
hydrodynamics of a ffuid made of ions and electrons with 


heat conduction and CRs diffusion is 

dp 


dt 

dpu 

de 

'dt 


V.{pu) =0, 

V. (^puu+ptot - 

V. ( (e -l-ptot)M - 


BB 


47r 

B(B.u) 


= 0 , 


47r 


= -V.F, 


cond 


V.F< 


CR ; 


dB 

dcE 

dt 

deer 

dt 


— X {u X B) = 0, 

+ V.(eEM) = —PeV.M — V.Fcond + 'HeI 
+ V.(ecrM) = -PcrV.M - V.FcR , 


( 1 ) 

( 2 ) 

(3) 

(4) 

( 5 ) 

( 6 ) 


where p is the gas mass density, u the gas velocity, B the 
magnetic field, e = 0.5pu^ + eth + Scr + is the to¬ 

tal energy density, eth = e\ + ce is the thermal energy 
density, ce and ej are the electron and ion energy densi¬ 
ties, Ccr is the energy density of CRs, and total pressure 
Ptot = (j - l)eth + (7cr - l)ecr + 0.5R^/47r, where 7 and 
7 cr are the adiabatic indexes of the gas and CRs, respec¬ 
tively. All energy components are energies per unit vol¬ 
ume Ci = Ei/Ax^, where Aa: is the cell size. The heat flux 
Tcond is carried by electrons alone, and the temperature 
of ions and electrons couple through the energy exchange 
term Hei- Since ions are heavier than electrons, their ther¬ 
mal velocity is lower by a factor oc •^/we/ttii- Therefore, 
for a given increase in thermal velocity at a shock, it cor¬ 
responds to a higher temperature increase for ions than for 
electrons. It justifies the fact that electron energy is treated 
in a separate energy equation that does not fulfil the jump 
condition. The same arguments apply to CRs compared to 
the thermal gas since velocities of CRs are much higher than 
thermal velocities of the ions (a.k.a. CRs are a non-thermal 
gas component). 

We assume that the thermal components of the gas are 
made only of ions and electrons, meaning that the gas is 
fully ionised. In some situations, this approximation will 
break down (star-forming regions), but it is a good approx¬ 
imation for the warm or hot phases of the ISM and for 
the intergalactic gas. CRs also diffuse through the diffusion 
flux term i^cR- Equation ([S]) can be repeated for as many 
CR energy components required to decompose the CR en¬ 
ergy power spectrum (with different diffusion coefficients). 
Introducing multiple CR energy bands requires including 
adiabatic and radiative losses for a proper treatment of the 
CR power spectrum, so we defer this to future work and 
simply assume that CRs are represented by one ffuid with 
one effective diffusion coefficient and one effective adiabatic 
index. 

W e use the AMR RAMSES code detailed in iTevssieil 
(|2002fl . Equations ([51) and ([51) for electrons and CRs, re¬ 
spectively, are added to the set of MHD equations. The full 
set of equations is solve d with the standard MHD solver 
of RAMSES described in iFromang et al.l (|200fff) , where the 
right-hand terms of equation m are treated separately as 
source terms. The ion internal energy is obtained for free by 
subtracting the electron internal energy from the total ther¬ 
mal energy. The induction equation (equation [p is solved 
using constrained transport (iTevssier et al.l2006f) . Godunov 
fluxes are modified to account for the extra energy compo- 
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nents and total pressure made of ions, electrons, and CRs. 
Therefore, the effective sound speed used for the CFL time 
step condition accounts for the extra pressure components 
(i.e. total pressure of the fluid). 

2.2. Anisotropic diffusion 

We explain the method for the conduction of heat through 
electrons, but this can be applied to a single tempera¬ 
ture model (ions and electrons at thermal equilibrium) or 
for the diffusion of CRs. The difference for CRs is that 
the diffusion coefficient Dck is constant throughout the 
simulation volume and that it applies to the gradient of 
CR energy instead of the gradient of temperature, e.g. 
Fen = —Dcnb{b.V)ecR, where b is the magnetic field unit 
vector. 

In the presence of magnetic fields, the conduction of 
heat in a plasma is 


den 

dt 


— V.Fcond = —V. ( —K||b(6.V)TE)—V. ( —KisoVRe) , 

(7) 


where Kiso and Ky are the isotropic conduction and parallel 
along the magnetic field lines coefficients respectively, with 
K|| = Ksp - Kiso, and Kiso = /isoKsp with /iso an arbitrary 
small coefficient for anisotropic diffusion (/iso = 1 is for 
isotropic diffusion). In many astrophysical cases, Kiso/Ky <C 
1 since the Larmor radius is much smaller than the mean 
free path of electrons. For instance, in the hot gas of galaxy 
clusters with temperature Te = 3 keV, electron density 
riE = 10“^cm“^ and B = 1 pG, the Larmor radius is 
Al = 10® cm and the mean free path Amfp = 10^^ cm. How¬ 
ever, to ensure numerical stability, we set up the isotropic 
conduction term to be around one per cent of the parallel 

conduction coefficient. _ 

Th e conduction coefficient for electrons is the ISoitzeil 
(|195(1^ value 


Ksp = nnknDc , (8) 

where kn is the Boltzmann constant and Dq the thermal 
diffusivity 


Dc = 8x 10®^ 


Te \ 
lOkeV ) 


He 


cuff s ^. 


(9) 


2.2.1. Implicit scheme 

This diffusion step is performed after the magnetohydrody¬ 
namics step. Here, we explain how the anisotropic conduc¬ 
tion and diffusion equation is solved for a uniform grid, and 
the AMR part of the solver is explained in section 12.2.21 
Given the restrictive stability condition for an explicit 
scheme on the equation dZD, Afdiff = Aa;^/(2R>), compared 
to the behaviour of the CFL condition for the hydrodynam¬ 
ics time step Ath = A.x/{u -f Cs), using an implicit solver 
on the equation diffusion will alleviate this numerical sta¬ 
bility constraint. Thus, the diffusion term is solved over one 
hydrodynamical time step using implicit integration. Dis- 
cretising equation © in ID leads to (we omit subscripts E 


and cond for clarity) 


e"+i 


At 


_ ^n+1 

2+4 i—i 


Ax 


= e,- 


( 10 ) 


where the subscript i is for the cell position, and the super¬ 
script n for the time. Each of the fluxes is evaluated 

at the cell interface (2 in ID, 4 in 2D, and 6 in 3D) at time 
n -|- 1. This can be rewritten as 


n+1 


-At 


n (rpn+l 

Gz+1 


^ 2 


Ax^ 


7;-V) 



where = riiksT /(q—1) {ui here is the gas number density, 
not time index n). Eg nation 1111 implicitly assumes that the 
cells are cubic, which is the case in RAMSES. We use the 
value of K at time n since we employ an implicit solver, 
which requires the constancy of the coefficient during the 
integration step. Thus, defining = K^^iAt/Ax'^ and 

Cv,i = nikn/i'^ — 1), we obtain 




n+1 /^n rpn+1 _ /^n 


( 12 ) 


We discretise equation © in 2D (3D can be obtained 
from 2D with little effort) assuming the cells have the same 
extent in x, y, z directions 


p"+i 


f: 


n+1 


At- 




^n+l 


^n+1 

i- 1 d 


^n+1 

id-k 


Ax 


^id 


(13) 


for cell position //. These quantities are e valuated with 
the ce ntred symmetric scheme proposed by iGiinter et al.l 
(I 2 OO 5 II for the anisotropic part of the flux. The anisotropic 
flux at cell interfaces F.’Rl . and are evaluated 

from their cell corner fluxes +“ 1/2 j±i/ 2 J 


TT'ani 

i+id ~ 


TT'ani 

^id+h ~ 


TT'ani 

T--r 2'J 2 


TT'ani 

^+h3+l 


TT'ani 


TT'ani 


The anisotropic cell corner flux is 


TT'ani 

^i+hd+l 


= Ki\b. 


II Ox 



dx dy) 


(14) 


where barred quantities are arithmetic averages over the 
cells connected to the corner; i.e., 


dT 

dx 

dT 

dy 


2 

+ 5 V,i-+l,3 + ^ 

2 


rpn-\-\ I rrin+l rpn-\-\ rrnn+l 

2 Ax 


T-in+l I T-in+1 T-in+1 rrin+l 

2Ax 


\\id ^ \\i+id ^ hd+i ^ ^Ib+I J+I 

4 
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For the isotropic part of the flux, we use a classical dis¬ 
cretisation, where the fluxes of equation m are simply 
obtained from the left and right states of the cell interface; 
i.e., 

riiso _ ^isoz+l,j ^i,j /-, 

^ • (15) 

With this numerical scheme, the matrix system Ax = c 
formed by equation (na, with A the matrix including the 
C" and C" coefRcients, x the vector of temperature 
and c the vector of energy densities can be gener¬ 

alised in multi-dimensions by equation m- The sparse 
matrix A is positive definite and symmetric, so we can use 
the conjugate gradient a lgorithm to solve t his sy stem of 
linearised equations as in ICommercon et al.l (|201lli (for ra¬ 
diation hydrodynamics in their case). 


2.2.2. Diffusion with adaptive mesh refinement 

We fo llow the procedure introduced in ICommercon et ahl 
^|2014^ for solving the diffusion equation with AMR and 
adaptive time-stepping on a level-by-level basis. Each level 
i is evolved with a time step that is twice less than 
the coarser level i — 1, At^ = At^“^/2, meaning that level 
i evolves with two consecutive time steps before doing one 
time step of level i — 1. The adaptive time-stepping is per¬ 
formed for all solvers in RAMSES including the diffusion 
solver presented in this paper, and finer levels are updated 
first. 

At a given level of refinement i, there are two types of 
non-uniform interfaces: the fine-to-coarse interface (inter¬ 
face between a cell at level i and a cell at i — 1) and the 
coarse-to-fine interface (interface between a cell at level £ 
and a cell at £ -|- 1). In both cases, we use Dirichlet bound¬ 
ary conditions: cell values at level boundaries are imposed. 
One could choose Neumann boundary conditions, which im¬ 
pose the fluxes and guarantee energy conservation as what 
is done for the hydrodynam ical solver. However, as shown 
in ICommercon et al.l (|2014ll . this sort of imposed flux con¬ 
ditions could lead to negative values for lang time steps. 

For the fine-to-coarse interface, we use values of £ — 1 
at time n as imposed boundary cond itions for level £. With 
the minmod scheme (|van Leeiill979|l (monotonised central 
is also a valid choice), we interpolate the values of level £—1 
on a finer virtual grid at level £ to impose the fine-to-coarse 
boundary. For the coarse-to-fine interface, we use values of 
£ -|-1 at time n -|-1 as imposed boundary conditions for level 
£. We restrict the value of the boundary coarse cell at level £ 
to the average value of the 2'^™ cells of level £ -I- 1 to impose 
the coarse-to-fine boundary. In any case, eq. [TS] remains 
correct since the level-by-level diffusion solver only deals 
with data estimated at the same level of refinement. The 
combination of the use of Dirichlet boundary conditions at 
level interfaces and of interpolation or restriction operations 
does not break the symmetry of matrix A. The imposed 
values of the neighbouring cells at different refinement levels 
go into the right-hand side (vector c) of the matrix system 
Ax = c, because it is the case for the computational domain 
imposed boundary conditions. 


2.2.3. Limitations of the method 


let boundary at level interfaces and the non-monotonicity 
preserving nature of the solver. To circumvent the first lim¬ 
itation, one could adopt a unique time step strategy, where 
diffusion is solved for all levels at once. This guarantees en¬ 
ergy conservation but slows down the calculation. A good 
compromise would be to do the diffusion step rather than 
every fine time step of the simulation, but only over coarse 
time steps (or even ever y 10, 100, etc, coarse tim e steps). 
Nevertheless, as shown in ICommercon et al.l (|2014f) . Dirich¬ 
let boundary conditions are robust in most cases, and in 
practice, energy conservation is acceptable (see tests be¬ 
low). 

Concerning the non-monotonicity preserving nature of 
the diffusion solver, since the flux at a cell interface (equa¬ 
tion I14|) includes a transverse gradient of temperature, it 
is not guaranteed that the flux flows from the high tem¬ 
perature cell to the low temperature cell. Thus, the tem¬ 
perature can become negative, and the monotonicity of the 
solution is not preserved (see for instance the test prob ¬ 
lem illustrated in figure 5 of ISharma fc Hammett! [2007f) . 
ISharma fc Hammetfi (l2007f) propose a slope limiter for the 
centred symmetric scheme with explicit time integration 
that preserves the monotonicity of the solution. However, 
employing a slope limiter would not allow us to use the 
conjugate gradient algorithm any longer since matrix A be¬ 
comes non-symmetric. In our case, we prefer to conserve the 
simple and fast solution offered by the conjugate gradient 
algorithm and to fix the negative temperatures created by 
the anisotropic conduction solver with tiny positive values. 
Future developments will include slope limiters, and the 
matrix system will be s olved using a bi-conju gate gradient 
algorithm instead (as in iGonzalez et aLll2015f) . 


2.3. Ion and electron temperature coupling 


We now explain how we model the coupling of the tem¬ 
perature of ions and electrons in RAMSES. This coupling 
step is performed after the diffusion step. The energy ex¬ 
change rate 'Hei introduced in equation ([5]) is written for 
a gas that we assume to be fully ionised and composed of 
hydrogen and helium 


dcE 

dt 

dei 

dt 


= T-Lei = 


= —Hei = 


Ti — Te re^b 
"Teq.EI 7 — 1 ’ 
Te — Ti nifee 


■^eq.IE 7—1 
with the equilibration timescales 



3 

SmETOifeg / 

^Te ^ 

’^eq,EI 

8(27r)5niZj^(7|. In A ' 

^mE 


3 

SmETOl^E / 

^Te\ 

and 

8(27r)5niZj^(7|. In A ' 

v^e/ 


3 

SmEmik^ 

/Te 

"^eqJE — 

8(27r) 5 riEZfq^ In A 

\mE 


CO 

(Te\ 


8(27r) 2 riEZfq^ In A 



IiV 

mi) 


]±\ 

mi) 


(16) 


(17) 


(18) 


There are two main limitations to our current anisotropic where mE and mi stand for the electron and ion mass, re¬ 
diffusion solver: the non-energy conservation due to Dirich- spectively; ni = p/{pimp) is the ion number density with 
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Trip the proton mass and being the ion mean molecular 
weight; Zi is the ion charge number; the electron charge; 
and In A = 40 is the Coulomb logarithm. Since ms <C mi, 
we neglect the Ti/mi term in the energy exchange rates. 
For a mixture of hydrogen and helium, the timescale of the 
temperature coupling between hydrogen and helium would 
be proportional to Thus, the equilibra¬ 

tion between hydrogen and helium operates 10 ® times faster 
than the equilibration timescale between hydrogen and elec¬ 
trons or helium and electrons. For hydrogen or helium, the 
ratio of mn.He/^H He simplihes to mp (not the case 
for heavier elements). Therefore, protons and Helll can be 
considered as one single population sharing the same tem¬ 
perature. 

Writing feq = Teq,Ei/?^E = Teq.iE/?^!, we See that 



0.0 0.2 0.4 0.6 0.8 1.0 


X 


deE _ dei _ Ti — Te fee 
dt dt feq 7 — 1 

The variation in energy is symmetric between the transfer 
of ion energy and electron energy, but this is not the case for 
the variation in temperature, which is dtT^j oc n^\dteE.i- 
Therefore, the difference in temperature change between 
ions and electrons will be a factor fii/fiE, where pE is the 
“mean molecular weight of electrons” fiE = p/{nEmp). For 
astrophysical plasmas essentially composed of hydrogen and 
helium, the change in temperature between ions and elec¬ 
trons is close to symmetric. 

To solve for the system of two non-linear coupled equa¬ 
tions, we rewrite HU) as 


Gi(r^+\rj"+i) 

G2(r^+\rj"+®) 


= At 




K[T^+^-TS] =0, 


= At 




- K [T”+i - Tj"] = 0 , 


( 20 ) 


where K = 3mpA:B®/[8(27rmE)° ®nE 9 E In A]. We introduce 
the vector X = (Tg'*'^, T"”''^)^, and rewrite the previous 
two equations as 

Gi(A'= + AA)-Gi(A'=) = AA, 

G2{X>^ + AX) - G2{X’^) = AA. 

The goal is now to iterate on AA to find the value of 

Gi(A^ -I- AA) = 0 and G 2 (A^ -|- AA) = 0 as required 

by equations (l20)) . The derivatives dxG\ and dxG\ are 
obtained by differentiating equations ( 1201 ) with respect to 
Tg and , where the superscript k replaces the superscript 
n -|-1 in (EQl). The values of AA are obtained with Cramer’s 
rule, and AA is added to the value of A^ until the variation 
in AA/A'^ < lO-'^. 


3. Numerical tests 

If not indicated, we use a courant factor of G = 0.8 with 
adaptive time-stepping between the different levels of re¬ 
finement. 


Fig. 1. Top: Diffusion of a pressure step function as a function 
of position at t = 9.3 x 10“^ (black), t = 1.9 x 10“® (red), t = 
5.6 X 10“® (blue) for the simulation (symbols) and the analytic 
solution (solid lines). The extra levels of refinement are indicated 
as dashed lines. Bottom: Relative errors of the simulation with 
respect to the analytic solution. 


3.1. Diffusion of a step function 

We performed a simple ID test of the diffusion of a step 
function with a constant diffusion coefficient of D = 1 
without hydrodynamics. The initial setup is a pressure of 
Pi = 0.8 between x = [0.25,0.75] and Pq = 0.4 outside 
with Gy = P^b/ (prup) = 1. We start with a minimum level 
of refinement of £min = 3 and refine up to level fmsLX = 7 
wherever the relative pressure variation is larger than 10 %. 
The time step is chosen to be At^ = 7 x Atuiff (for level 
7), where Atuiff is the stability time step for diffusion (only 
useful for explicit schemes). 

The time evolution of the step function has the following 
analytical solution for a constant diffusion coefficient 

P(x. t) = , (21) 

where xq = 0.25. Figure [1] shows the result of the diffusion 
of the step function with the analytical solution at three 
different times. The analytical solution is reproduced well, 
below a 10% relative error. 

3.2. Temperature coupling 

Figured] shows the result of the one cell problem with two 
different initial temperatures for electrons Tep = 10^° K 
and ions Tj o = 10® K. The solid lines indicate the solution 
with a time step of Teq_Ei,o/50i ^^nd points are the solution 
with a time step of 20Teq,Ei,o- They are in excellent agree¬ 
ment, even for the run where the time step is 20 times the 
timescale of temperature coupling. 


3.3. Sod test with two energy components 

We set up the conditions for a ISodI (Il978ll shock tube with 
two temperatures and without conduction or temperature 
coupling. The initial left and right gas density is pn = 1, 
Pr = 0.125, the first pressure Piy = 0.34, P 13 = 0.066 (gas 
or ion), the second pressure P 2 ,l = 0.66, P 2 ,r = 0.034 (re¬ 
spectively CR or electron), thus Pl = 1 and Pr = 0.1, and 
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log time (Myr) 

Fig. 2. Electron temperature (red or diamonds) and ion tem¬ 
perature (blue or pluses) as a function of time for the one cell 
problem. The solution for a time step of req,Ei,o/50 is represented 
in solid lines and for the time step of 20req,Ei,o with symbols. 






Fig. 3. Solution of the lSodI lll978l l shock tube problem with two 
temperatures. There is no diffusion and coupling of the temper¬ 
atures. We plot the velocity (top left), pressure (bottom left), 
and density (bottom right) at t = 0.245. The squares correspond 
to the result of the simulation and the red lines to the analytical 
solution. The light blue and dark blue squares in the pressure 
plot stand for the first and second species, respectively, while the 
black squares are for the total pressure. The level of refinement 
is plotted in the top right panel. 


zero velocity. The adiabatic index for the first and second 
temperatures is 71 2 = 1.4. The minimum level of refine¬ 
ment is ^min = 3, and we refine up to level imax = 10 wher¬ 
ever the relative gradient of any hydro variable is greater 
than 5%. Figure [3] shows the result of the simulation, which 
exhibits extremely good agreement with the analytical so¬ 
lution. At the contact discontinuity x = 0.7, the total pres¬ 
sure is continuous, but not the individual pressure terms. 
We repr oduce the structure of a two-component shock well 
(see e.g. iPfrommer et al.ll200(il for comparison). 

3.4. Shock with two temperatures and heat conduction 

We model the formation of a ID stationary shock typically 
arising at the virial radius of a cluster of galaxies. The gas is 
of primordial cosmological composition with a 0.76 fraction 
of hydrogen and 0.24 of helium and is fully ionised. Thus, 


the mean molecular weight of ions, electrons, and gas are 
Pi = 1.22, pe = 1.13, and p = 0.59. The right state is 
with gas density ur = 10 “'^ cm“^, gas velocity ur = 2000 
km s“^ oriented in the negative x-direction, and ion and 
electron temperatures at equilibrium Ti,r = Te,r = Tr = 
4 X 10^ K. Using the Rankine-Hugoniot jump conditions, 
we find that the corresponding left state is np — 2.3 x 
10“^ cm“^, ul = 875 km s“^, and Ti,r = Te^r = Tr ~ 
8.4 X 10^ K for an adiabatic index of 7 = 5/3. The box 
length is 8 Mpc, with a minimum and a maximum level 
of refinement of fmin = 4 and £max = 9, respectively, and 
with refinement triggered wherever the relative gradient of 
any hydro variable is larger than 10%. The left and right 
boundaries are imposed with the initial values of left and 
right states. 

Figured shows the result of the simulation with the two 
temperatures including conduction of the electrons and the 
temperature coupling between ions and electrons, as well 
as the solution for a single temperature with or without 
conduction. For the single temperature experiment without 
conduction, the initial left and right states are preserved at 
time t = 10 Gyr. With the addition of conduction, a ther¬ 
mal precursor forms ~ 1 Mpc ahead of the shock, which is 
present for both the single and the two-temperature exper¬ 
iments. 

In the case of the two temperatures, we see that the 
electrons have the higher temperature value in the ther¬ 
mal precursor compared to the ions with a slow rise similar 
to the thermal precursor in the single-temperature experi¬ 
ment. Ion temperature in the precursor lies below—with a 
sharp rise at the shock interface—because the coupling of 
the two temperatures is not instantaneous. The increase in 
ion temperature right after the shock interface compared to 
the average gas temperature is because the ions absorb the 
upstream kinetic energy at the shock, and electrons only 
lag behind and thermally recouple with the ions. Far away 
from the shock front, bo th temperatures are at equi librium. 
The reader can refer to IZel’dovich fc Raized (^1967^ for the 
detailed description of a two-temperature shock with con¬ 
duction. 


3.5. Conduction in a circular loop 

Anisotropy can be tested by the propagation of heat 
thr ough a magnetic loop, here in 2D (see for compariso n 
e.g. ISharma fc Hammett) 1200^ iRasera fc Chandranll200§) . 

We initialised a patch of temperature of T = 10® K wher¬ 
ever 2 < r < 3 kpc and 165 < 9 < 195° in a background 
gas with temperature T = 10® K and density n = 1 cm“®. 
The size of the box is 10x10 kpc, and we performed two 
runs either with a uniform grid 128^ or with AMR with 
•^min = 4 and ^max = 8 and refinement triggered if the rel¬ 
ative gradient of pressure is larger than 20%. We used a 
Spitzer conductivity with K || = 0.99/tSp and kiso = O.OlKsp- 

The result of this numerical experiment is shown in 
Fig. [S] for two runs with and without AMR. The initial 
patch of temperature has diffused anisotropically along the 
circular-loop shape of the magnetic field lines. The results 
of the AMR and of the uniform grid runs are in very good 
agreement, as shown by the white and red contours. 

We also explored the effect of varying the perpendicu¬ 
lar conductivity coefficient with Kiso = [0.001, 0.01, 0.1]/csp- 
The resulting temperature maps are shown in Fig. As 
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Fig. 4. Stationary cosmological shock experience including conduction and two-temperature coupling. The result of the two- 
temperature run (blue) is represented with dotted lines for the ion temperature (pressure), in drished for the electron temperature 
(pressure), and solid for the average temperature (total pressure) of the gas at f = 10 Gyr. For comparison, the results of the 
single temperature runs are plotted in red with conduction and in black without conduction. An electron precursor forms in the 
upstream part of the shock front. 



2 kpc 2 kpc 


Fig. 5. Temperature map at t = 2 Gyr for the heat conduction 
in a 2D magnetic loop without AMR (grey levels). The solution 
with AMR is plotted with red contours for T = 5 x 10® K and 
T = 10^ K (white contours without AMR). 


expected, the spread of temperature perpendicular to the 
magnetic field lines increases with the increase in the 
isotropic component of the diffusion coefficient. 

3.6. Sovinec test 

ISovinec et al.l (I2004D designed a test to measure the numer¬ 
ical perpendicular diffusion Knum for the anisotropic diffu- 


Fig. 6. Temperature map at t = 2 Gyr for the heat conduc¬ 
tion in a 2D magnetic loop with AMR and different fractions 
of perpendicular conductivity coefficients kiso = 0.1 KSp (red 
contours), Kiso = O.Olftsp (grey levels and white contours), 
and Kiso = O.OOIksp (blues contours). Two contours are plot¬ 
ted for each simulation result corresponding to temperatures 
T = 5 X 10® K and T — 10^ K (except for kiso = O.Iksp where 
only T = 5 X 10® K is reached). 

sion in 2D. The energy evolution follows 
de 

— = -V. (-K||b.(b.V)rE) - V. (-Ki.oVTE) + Q , (22) 

where Q is a heating source term of the form Q(x, y) = 
Qo cos(7ra;) cos(7ry). The magnetic field is chosen such that 
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A X 

Fig. 7. Results of the perpendicular numerical d iffusion coeffi¬ 
cient fi:num as a function of resolution Aa: for the lSovinec et ahl 
(l2004fl test problem in a solid line with square symbols. The 
numerical diffusion scales with Ax^ 

(dotted line). 

B.VTe = 0, i.e. the parallel heat flux component is 
suppressed, and only the perpendicular component re¬ 
mains. The analytical stationary solution at the centre is 
T(0,0) = Qo/[27r2(Kiso + Knum)]- For this test kiso = 1 
(kSp = 100), p = 1, Qo = 27r2, = cos(7ra:) sin(7ry), 

and By = — sinfTra:) cos(7rp). The box size ranges ove r 
[—0.5,0.5] X [—0.5,0.5]. As in lR.asera & ChandranI (|2008ri . 
we imposed the temperature at the boundaries of the sim¬ 
ulation box with Tbound = 0. 

We ran the simulation for various uniform grid res¬ 
olutions from 16^ to 256^ up to the steady-state solu¬ 
tion, and we measured th e valu e of rmea(0,0). We fol¬ 
lowed ISharma fc HammettI (|2007ll to obtain the value of 
Knum as follows K„um = Tiso(0,0)/Tmea(0,0) - 1, where 
riso(0,0) is the the central temperature calculated in the 
isotropic limit {kiso = ^Sp)- This measurement is more 
accurate than assuming that the isotropic diffusion gives 
riso(0,0) = 1 exactly. The result is shown in Fig. [71 and we 
see that the perpendicular numerical diffusion coefficient 
scales quadratically with resolution oc Ax^ and that, even 
at low resolution Ax = 1/16, the numerical diffusion is a 
factor 100 below our explicit perpendicular diffusion coeffi¬ 
cient Kiso- 

3.7. Supernova explosion in a 3-phase medium: ions, 
electrons, and cosmic rays 

To fully validate the implementation of anisotropic diffusion 
of CRs and conduction of heat through electrons and their 
coupling with ions in an astrophysical context, we model 
the explosion of a supernova in 3D. The explosion is trig¬ 
gered in a homogeneous medium of density n = 1 cm“^ 
and with temperatures of electrons and ions at equilibrium 
T = 10^ K, and the background energy density of CRs is 
one-third of the total energy. We assume that the gas is 
totally ionised and composed of a cosmic mixture of hydro¬ 
gen and helium, thus, pi = 1.22, /is = 1.13, and p = 0.59. 
The magnetic field is uniform along the x-axis of the box 
and with a magnitude of 1 pG, thus the ratio of internal to 
magnetic pressure is initially of /3 = 45. We inject in the 
eight central cells a total of Rsn = 10^^ erg of energy with 


1/3 in CR energy, 1/3 in ion energy, and 1/3 in electron 
energy. 

The temperature of electrons is conducted at the Spitzer 
rate and couples with the ion temperature. Both ion and 
electron adiabatic index are equal to 7 = 5/3. CRs are 
diffused with a diffusion coefficient Dcr = 3 x 10^^ cm^ s“^ 
and have an adiabatic index of ycR =4/3 corresponding to 
nltra-relativistic particles. As cnrrently implemented, our 
solver could allow for several CR energy components with 
various diffusion coefficients and CR adiabatic indexes. We 
nse a box size of 500 pc with a coarse grid of level .^min = 5 
and refine up to £„iax = 9 wherever the relative gradient of 
any hydro variable is more than 20 %. 

For the first test, we did electron and ion coupling with- 
ont diffusion or conduction, and the second test was with 
isotropic diffusion and conduction. Figure [5] shows the result 
of the two simulations after time f = 4 kyr and is compared 
to the analytic prediction from I Sedov! (|1959f) . In the ab¬ 
sence of conduction and diffusion, the analytical solution is 
reproduced well, except that the shock interface is smeared 
by the numerical diffusion. The electrons and ions are not 
yet at thermal equilibrinm within the bnbble. The solution 
is slightly modified with isotropic diffusion of CR energy 
and conduction of electron temperature. The gas density 
exhibits two peaks: one ahead of the shock front position 
and another behind. The ahead over-density is the product 
of the conduction of electrons over a few parsec, while the 
second over-density lags behind. The total pressure within 
the bubble is decreased since the CR and electron energies 
have leaked out of the bubble. Figure [H] shows the result at 
time t = 100 kyr. Again, the experiment withont conduc¬ 
tion and diffusion reproduces the analytical prediction well, 
and ions and electrons are at thermal equilibrium. At this 
time, the Sedov shock has taken over the electron precur¬ 
sor, which has disappeared. It turns out that electrons and 
ions are at thermal eqnilibrium. However, a large amonnt 
of the CR energy has escaped the bubble, so that there is 
less pressure within the bubble, explaining why the shock 
front slightly lags behind. 

The final experiment was run with anisotropic conduc¬ 
tion of electrons and anisotropic diffusion of CRs. FigurefTUl 
shows 2D slices in the middle xz plane of several quanti¬ 
ties: temperatures, density, CR energy, velocity, and mag¬ 
netic field amplitude at two different times t = 5 kyr and 
t = 100 kyr. We see that the electron temperature and CR 
energy are diffused along the a;-axis. The ion temperature 
catches up with the electron temperature in the heat pre¬ 
cursor, but it takes some time for the two species to reach 
thermal equilibrium. A jet-like nozzle forms in the horizon¬ 
tal direction due to the anisotropic over-pressurised region, 
although it propagates more slowly than the genuine shock 
front. Nevertheless, at later times, the explosion approaches 
spherical with still some degree of anisotropy in the bub¬ 
ble expansion. In particular, the shock front moves faster 
along magnetic field lines (x-axis) than perpendicular to 
them (note the amplification of the magnetic field due to 
the compression along the 2 ;-axis). This results in an aspher- 
ical bnbble shape, which is in advance in the x-direction but 
behind in the ^-direction (as for the y-direction). 

4. Conclusion 

We have presented the implementation of a new solver 
for the anisotropic diffusion along magnetic field lines of 
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Fig. 8. Spherically averaged values of density (left), pressure (middle), and velocity (right) as a function of radius at t = 4 kyr 
for the 3D Sedov explosion without diffusion (top), and with isotropic diffusion (bottom). The solid lines are the result of the 
simulation, and the dashed line is the analytic prediction for the standard Sedov solution. The pressure is decomposed into total 
pressure (black), ion pressure (blue), electron pressure (red), and CR pressure (magenta). Without diffusion, the analytic solution 
is reproduced well. With diffusion, there are both a CR precursor and an electron precursor. 




r (pc) r (pc) 



r (pc) 





r (pc) r (pc) r (pc) 


Fig. 9. Same as Fig. [H] at t = 100 kyr. The shock front takes over the conduction front of electrons. A CR pressure precursor is 
still visible in the conduction and diffusion case with a subtle effect on the Sedov shock front position, which slightly lags behind 
the solution without conduction and diffusion. 
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Fig. 10. Map cuts in the xz plane for the 3D Sedov explosion including anisotropic electron heat conduction and anisotropic CR 
diffusion at t = 5 kyr (first two rows) and t — 100 kyr (last two rows). From left to right and top to bottom: the ratio of electron 
to ion temperature, the ion temperature, the electron temperature, the gas density, the magnetic field amplitude, the gas velocity 
amplitude, the ratio of CR to internal energy, and the CR energy density. The black dotted circles are the shock front positions of 
the standard Sedov solution R{t) ~ 1.1527(FsN/p)^^^t^^®. The explosion exhibits a nozzle-like heat/CR precursor in the horizontal 
direction due to the anisotropic nature of the diffusion (along the horizontal magnetic field). 


CRs and of heat with two-temperature components, elec¬ 
trons, and ions in the AMR code RAMSES. This solver is 
implicit in time and supports adaptive-time-stepping on 
an AMR level-by-level basis. The implicit anisotropic dif- 
fusio n scheme uses a conjugate grad ient method adapted 
from ICommercon et ahl (|201lL 1201411 . We tested this new 
implementation against several basic numerical experi¬ 
ments with, for some cases, typical astrophysical conditions, 


as well as a first simple application to a supernova explosion 
in a uniform medium. This preliminary work indicates that 
anisotropic conduction or diffusion, together with tempera¬ 
ture coupling, is relevant to several astrophysical problems, 
as already highlighted in the literature, and, in particular, 
may affect the propagation of remnants of supernovae. 

Future applications of this new part of the RAMSES code 
will include the study of active galactic nuclei jets in the 
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intra-cluster medium, virial shocks in cosmological simu¬ 
lations of galaxy clusters, CR-driven galactic winds, su¬ 
pernovae explosions, and thermal instabilities in the ISM, 
amongst others. 
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